/*
 * HelloItpp_test.cpp
 *
 *  Created on: Jan 5, 2013
 *      Author: shenojia
 */

#include "HelloItpp.h"

#include <stdio.h>
#include <stdarg.h>
#include "CppUTest/TestHarness.h"
#include "cppunit/TestAssert.h"

TEST_GROUP(HelloItpp)
{
	void setup()
	{
	}
	void teardown()
	{
	}
};

IGNORE_TEST(HelloItpp, Hello)
{
	printf("hello jiagui, welcome back to itpp exercise\n");
}

IGNORE_TEST(HelloItpp, Vec)
{
  /*-----------------------------------------------------------------------------
   *  test of vector
   *-----------------------------------------------------------------------------*/
  cout << "=============================" << endl
       << "   Testing Vec<bin> (bvec)" << endl
       << "=============================" << endl;
  bvec bv1 = randb(10);
  bvec bv2 = randb(10);
  bin bx = randb();
  vec_common_operators(bv1, bv2, bx);
  vec_logical_operators(bv1, bv2, bx);

  cout << "=============================" << endl
       << "   Testing Vec<int> (ivec)" << endl
       << "=============================" << endl;
  ivec iv1 = randi(10, 1, 9);
  ivec iv2 = randi(10, 1, 9);
  int ix = randi(1, 9);
  vec_common_operators(iv1, iv2, ix);
  vec_logical_operators(iv1, iv2, ix);

  cout << "===============================" << endl
       << "   Testing Vec<double> (vec)" << endl
       << "===============================" << endl;
  vec dv1 = randu(10);
  vec dv2 = randu(10);
  double dx = randu();
  vec_common_operators(dv1, dv2, dx);
  vec_logical_operators(dv1, dv2, dx);

  cout << "===============================================" << endl
       << "   Testing Vec<std::complex<double> > (cvec)" << endl
       << "===============================================" << endl;
  cvec cv1 = randn_c(10);
  cvec cv2 = randn_c(10);
  complex<double> cx = randn_c();
  vec_common_operators(cv1, cv2, cx);
  cout << "(a == c) = " << (cv1 == cx) << endl;
  cout << "(a != c) = " << (cv1 != cx) << endl << endl;


  // Test vectror initialisation with string
  vec v = "23.3 1232.7 0.111 1.525 0.333";
  cout << "Testing double vector initialisation with: "
    "\"23.3 1232.7 0.111 1.525 0.333\":" << endl << "v = " << v << endl;

  v = "-10.000 :.5:-4.5  1.33e+1, -.9, 1e0:1.5:1E+1";
  cout << "Testing double vector initialisation with: "
    "\"-10.000 :.5:-4.5  1.33e+1, -.9, 1e0:1.5:1E+1\":" << endl << "v = "
       << v << endl;

  ivec iv = "0xA :-0x1: -010";
  cout << "Testing int vector initialisation with: \"0xA :-0x1: -010\":"
       << endl << "iv = " << iv << endl;

  iv = "-5:3:9, 7, 1:10";
  cout << "Testing int vector initialisation with: \"-5:3:9, 7, 1:10\":"
       << endl << "iv = " << iv << endl;

  svec sv = "3 0xF -10, 0133 0177, 0x0 ";
  cout << "Testing short int vector initialisation with: \"3 0xF -10, 0133 0177, 0x0 \":"
       << endl << "sv = " << sv << endl;

  cvec cv = " (0.3, 0.4)  .2-.01i, 1e-3+0.25i";
  cout << "Testing complex vector initialisation with: \" (0.3, 0.4)  .2-.01i, 1e-3+0.25i\":"
       << endl << "cv = " << cv << endl;

  bvec bv = "1 1 0,1  1  ,  0 ,1  ";
  cout << "Testing bit vector initialisation with: \"1 1 0,1  1  ,  0 ,1  \":"
       << endl << "bv = " << bv << endl << endl;

  // Test of rem:
  v = "1.0 2.0 3.4 -4.5 6.7";
  double y = 0.76;
  cout << "v = " << v << endl;
  cout << "y = " << y << endl;
  cout << "rem(v, y) = " << rem(v, y) << endl;
  cout << "rem(10, v) = " << rem(10, v) << endl;
  mat M = "1.0 2.3; 4.5 -6.7";
  cout << "M = " << M << endl;
  cout << "rem(M, y) = " << rem(M, y) << endl;
  cout << "rem(10, M) = " << rem(10, M) << endl << endl;

  // Test of all and any:
  bvec b1 = "0 0 0 0 0 0 0 1 0 0";
  bvec b2 = "0 0 0 0 0 0 0 0 0 0";
  bvec b3 = "1 1 1 1 1 1 1 1 1 1 1 1 1";
  bvec b4 = "1 1 1 1 1 1 1 1 1 1 1 0 1";

  cout << "any(b1) = " << any(b1) << endl;
  cout << "any(b2) = " << any(b2) << endl;
  cout << "all(b3) = " << all(b3) << endl;
  cout << "all(b4) = " << all(b4) << endl;


}

IGNORE_TEST(HelloItpp, Mat)
{

  /*-----------------------------------------------------------------------------
   *  test of matrix
   *-----------------------------------------------------------------------------*/
  cout << "=============================" << endl
       << "   Testing Mat<bin> (bmat)" << endl
       << "=============================" << endl;
  bmat bM1 = randb(3, 4);
  bmat bM2 = randb(3, 4);
  bmat bM3 = randb(4, 3);
  bvec bv1 = randb(3);
  bvec bv2 = randb(4);
  bin bx = randb();
  mat_common_operators(bM1, bM2, bM3, bv1, bv2, bx);

  cout << "=============================" << endl
       << "   Testing Mat<int> (imat)" << endl
       << "=============================" << endl;
  imat iM1 = randi(3, 4, 1, 9);
  imat iM2 = randi(3, 4, 1, 9);
  imat iM3 = randi(4, 3, 1, 9);
  ivec iv1 = randi(3, 1, 9);
  ivec iv2 = randi(4, 1, 9);
  int ix = randi(1, 9);
  mat_common_operators(iM1, iM2, iM3, iv1, iv2, ix);

  cout << "===============================" << endl
       << "   Testing Mat<double> (mat)" << endl
       << "===============================" << endl;
  mat dM1 = randn(3, 4);
  mat dM2 = randn(3, 4);
  mat dM3 = randn(4, 3);
  vec dv1 = randn(3);
  vec dv2 = randn(4);
  double dx = randn();
  mat_common_operators(dM1, dM2, dM3, dv1, dv2, dx);

  cout << "==========================================" << endl
       << "   Testing Mat<complex<double> > (cmat)" << endl
       << "==========================================" << endl;
  cmat cM1 = randn_c(3, 4);
  cmat cM2 = randn_c(3, 4);
  cmat cM3 = randn_c(4, 3);
  cvec cv1 = randn_c(3);
  cvec cv2 = randn_c(4);
  complex<double> cx = randn_c();
  mat_common_operators(cM1, cM2, cM3, cv1, cv2, cx);


  cout << "========================================" << endl;
  cout << "   Testing initialisation with string" << endl;
  cout << "========================================" << endl;

  cout << "bmat M = \" 1 1 0; 0 1; 1 1 ,1,   1; ; 0 1\"" << endl;
  bmat bM = " 1 1 0; 0 1; 1 1 ,1,   1; ; 0 1";
  cout << "M =\n" << bM << endl;

  cout << "smat M = \"0xFF, -021 ,   100; 0,-0x01; 0xA, 10 012;  \"" << endl;
  smat sM = "0xFF, -021 ,   100; 0,-0x01; 0xA, 10 012;  ";
  cout << "M =\n" << sM << endl;

  cout << "imat M = \"0xFAC0, -021, 100000; 0,-0x01; 0xA, 10 012; \"" << endl;
  imat iM = "0xFAC0, -021, 100000; 0,-0x01; 0xA, 10 012; ";
  cout << "M =\n" << iM << endl;

  cout << "mat M = \".77 1e9 35; 0x7 3.5 -1000 \"";
  mat M = ".77 1.89e5 35; 7 3.5 -1000 ";
  cout << "M =\n" << M << endl;

  cout << "cmat M = \" 1.5+3i, (.33,1) ;  (333,-1) 2-0.2E-3i\"" << endl;
  cmat cM = " 1.5+3i, (.33,1) ;  (333,-1) 2-0.2E-3i";
  cout << "M =\n" << cM << endl << endl;


}
TEST(HelloItpp, MatFunc)
{

  /*-----------------------------------------------------------------------------
   *  test of Matrix fuction
   *-----------------------------------------------------------------------------*/
  cout << "=================================" << endl;
  cout << "    Test of matfunc routines     " << endl;
  cout << "=================================" << endl;

  cout.setf(ios::fixed);
  cout.precision(3);

  vec a = randn(5);
  cout << "a = " << a << endl;
  cout << "sum(a) = " << sum(a) << endl;
  cout << "cumsum(a) = " << cumsum(a) << endl;
  cout << "prod(a) = " << prod(a) << endl;
  cout << "sum_sqr(a) = " << sum_sqr(a) << endl << endl;

  mat A = randn(5, 5);
  cout << "A = " << A << endl << endl;

  cout << "sum(A) = " << sum(A) << endl;
  cout << "sum(A,1) = " << sum(A, 1) << endl;
  cout << "sum(A,2) = " << sum(A, 2) << endl << endl;

  cout << "cumsum(A) = " << cumsum(A) << endl;
  cout << "cumsum(A,1) = " << cumsum(A, 1) << endl;
  cout << "cumsum(A,2) = " << cumsum(A, 2) << endl << endl;

  cout << "prod(A) = " << prod(A) << endl;
  cout << "prod(A,1) = " << prod(A, 1) << endl;
  cout << "prod(A,2) = " << prod(A, 2) << endl << endl;

  cout << "sum_sqr(A) = " << sum_sqr(A) << endl;
  cout << "sum_sqr(A,1) = " << sum_sqr(A, 1) << endl;
  cout << "sum_sqr(A,2) = " << sum_sqr(A, 2) << endl << endl;

  cout << "repmat(a, 1, 3) = " << repmat(a, 1, 3) << endl;
  cout << "repmat(a, 3, 1, true) = " << repmat(a, 3, 1, true) << endl;
  cout << "repmat(A, 2, 2) = " << repmat(A, 2, 2) << endl << endl;

  cout << "Kronecker test" << endl;
  mat X = to_mat(randi(2, 2, 1, 4));
  mat Y = randn(3, 3);
  cout << "X = " << X << endl;
  cout << "Y = " << Y << endl;
  cout << "kron(X, Y) = " << kron(X, Y) << endl << endl;

  cout << "sqrtm of a real matrix" << endl;
  A = randn(3, 3);
  cmat A_sqrtm = sqrtm(A);
  cout << "A = " << A << endl;
  cout << "norm(sqrtm(A) * sqrtm(A) - A) = "
       << round_to_zero(norm(A_sqrtm * A_sqrtm - to_cmat(A)), 1e-13)
       << endl << endl;

  cout << "sqrtm of a complex matrix" << endl;
  cmat B = randn_c(3, 3);
  cmat B_sqrtm = sqrtm(B);
  cout << "B = " << B << endl;
  cout << "norm(sqrtm(B) * sqrtm(B) - B) = "
       << round_to_zero(norm(B_sqrtm * B_sqrtm - B), 1e-13) << endl << endl;

  cout << "Rank test" << endl;
  A = randn(3, 3);
  cout << "A = " << A << endl;
  cout << "rank(A) = " << itpp::rank(A) << endl;
  A.set_row(1, 3.0 * A.get_row(0));
  cout << "A2 = " << A << endl;
  cout << "rank(A2) = " << itpp::rank(A) << endl;
  B = randn_c(3, 3);
  cout << "B = " << B << endl;
  cout << "rank(B) = " << itpp::rank(B) << endl;
  B.set_col(1, B.get_col(0));
  cout << "B2 = " << B << endl;
  cout << "rank(B2) = " << itpp::rank(B) << endl;


}

IGNORE_TEST(HelloItpp,Array)
{
	Array<int> A1(10), A2(15);

	A1 = 1;
	A1(2) = 42;
	A2 = 5;

	cout << "Testing simple Arrays of integers:" << endl
	<< "A1 = " << A1 << endl
	<< "A2 = " << A2 << endl << endl;

	// Test of Array initialisation by string:
	Array<bvec> A3 = "{[1 1] [1 0 1 0] [0 0 1]}";
	cout << "Testing Array initialisation with: \"{[1 1] [1 0 1 0] [0 0 1]}\":"
	<< endl << "A3 = " << A3 << endl << "A3.lenght()" << A3.length()<< endl;

	bvec va;
	bvec vb = A3.shift_left(va);
	// Test of Array left_shift
	cout << "Testing of Array left_shift"
	<< endl << "shift out" << vb
	<< endl << "A3 after shift" << A3 << endl;

	Array<Array<imat> > A4 = "{{[5 3; 020 4] [1 0; 3 9]} {[0 -3; 1 0xa]}}";
	cout << "Testing Array initialisation with: \"{{[5 3; 020 4] [1 0; 3 9]} {[0 -3; 1 0xa]}}\":"
	<< endl << "A4 = " << A4 << endl << endl;

	// Test of operator()
	cout << "Testing Array::operator():" << endl
	<< "A4(1) = " << A4(1) << endl
	<< "A4(0)(1) = " << A4(0)(1) << endl << "A4.lenght()" << A4.length() << endl;

	// Test of left(), right() and mid() methods:
	cout << "Testing Array::left(), Array::right() and Array::mid():" << endl
	<< "A1.left(4) = " << A1.left(4) << endl
	<< "A1.right(5) = " << A1.right(5) << endl
	<< "A1.mid(2, 3) = " << A1.mid(2, 3) << endl << endl;

	// Test of swap function
	A4.swap(0, 1);
	cout << "Testing A4.swap(0, 1):" << endl
	<< A4 << endl;
}

IGNORE_TEST(HelloItpp,Mutual)
{
 //general parameters
    vec sigmaA = "0.01 2 4 7";//standard deviation (sqrt(variance)) of the mutual a priori information
    double threshold_value = 50;
    string map_metric="maxlogMAP";
    ivec gen = "07 05";//octal form
    int constraint_length = 3;
    int nb_blocks_lim = 10;
    int perm_len = int(itpp::pow10(3.0));//total number of bits in a block (with tail)
    double EbN0_dB = 0.8;
    double R = 1.0/3.0;//coding rate of PCCC
    double Ec = 1.0;//coded bit energy

    //other parameters
    vec sigma2A = sqr(sigmaA);
    int sigma2A_len = sigma2A.length();
    int nb_bits = perm_len-(constraint_length-1);//number of bits in a block (without tail)
    double sigma2 = (0.5*Ec/R)*pow(inv_dB(EbN0_dB), -1.0);//N0/2
    double Lc = -2/sigma2;//normalisation factor for intrinsic information (take into account the BPSK modulation)
    bvec bits(nb_bits);
    bvec tail;
    bvec bits_tail(perm_len);
    bmat parity_bits;
    int coded_bits_len = 2*perm_len;
    bvec coded_bits(coded_bits_len);
    vec mod_bits(coded_bits_len);
    vec rec_sig(coded_bits_len);
    vec intrinsic_coded(coded_bits_len);
    vec intrinsic_coded_p(2*nb_bits);
    intrinsic_coded_p.zeros();
    vec apriori_data(perm_len);
    vec extrinsic_coded;
    vec extrinsic_data;
    vec apriori_mutual_info(sigma2A_len);
    vec extrinsic_mutual_info(sigma2A_len);
    vec extrinsic_mutual_info_p(sigma2A_len);
    extrinsic_mutual_info.zeros();
    extrinsic_mutual_info_p.zeros();
    register int en,n,nb_blocks;

    //Recursive Systematic Convolutional Code
    Rec_Syst_Conv_Code rsc;
    rsc.set_generator_polynomials(gen, constraint_length);//initial state should be the zero state

    //BPSK modulator
    BPSK bpsk;

    //AWGN channel
    AWGN_Channel channel;
    channel.set_noise(sigma2);

    //SISO module
    SISO siso;
    siso.set_generators(gen, constraint_length);
    siso.set_map_metric(map_metric);

    //EXIT chart
    EXIT exit;

    //Randomize generators
    //RNG_randomize();
    RNG_reset(12345);

    cout << "=============================================" << endl;
    cout << "           Starting Simulation               " << endl;
    cout << " Extrinsic Mutual Information (IE) as a function of A priori Mutual Information (IA) " << endl;
    cout << "=============================================" << endl;
    cout << "  Block length = " << perm_len << endl;
    cout << "  Generator polynomials = " << std::oct << gen << std::dec << endl;
    cout << "  Turbo encoder rate 1/3 (plus tail bits)" << endl;
    cout << "=============================================" << endl;

    //main loop
    for (en=0; en<sigma2A_len; en++)
    {
        apriori_mutual_info(en) = exit.apriori_mutual_info(sigma2A(en));//a priori mutual info
        for (nb_blocks=0; nb_blocks<nb_blocks_lim; nb_blocks++)
        {
            //bits generation
            bits = randb(nb_bits);

            //RSC code
            rsc.encode_tail(bits, tail, parity_bits);//tail is added

            //form coder output
            bits_tail = concat(bits, tail);
            for (n=0; n<perm_len; n++)
            {
                coded_bits(2*n) = bits_tail(n);//systematic output
                coded_bits(2*n+1) = parity_bits(n,0);//parity output
            }

            //BPSK modulation (1->-1,0->+1)
            mod_bits = bpsk.modulate_bits(coded_bits);

            //AWGN channel
            rec_sig = channel(mod_bits);

            //first SISO RSC module  (tail is added)
            //intrinsic info. of coded bits
            intrinsic_coded = Lc*rec_sig;

            //a priori info. of data bits
            apriori_data = exit.generate_apriori_info(bits_tail);

            //SISO RSC module
            siso.rsc(extrinsic_coded, extrinsic_data, intrinsic_coded, apriori_data, true);

            //threshold
            extrinsic_data = SISO::threshold(extrinsic_data, threshold_value);

            //extrinsic mutual info
            extrinsic_mutual_info(en) += exit.extrinsic_mutual_info(extrinsic_data.left(nb_bits), bits);

            //second SISO RSC module (no tail added)
            //intrinsic info. of coded bits
            for (n=0; n<nb_bits; n++)
                intrinsic_coded_p(2*n+1) = Lc*rec_sig(2*n+1);//parity bits only

            //a priori info. of data bits
            apriori_data = exit.generate_apriori_info(bits);

            //SISO RSC module
            siso.rsc(extrinsic_coded, extrinsic_data, intrinsic_coded_p, apriori_data, false);

            //threshold
            extrinsic_data = SISO::threshold(extrinsic_data, threshold_value);

            //extrinsic mutual info
            extrinsic_mutual_info_p(en) += exit.extrinsic_mutual_info(extrinsic_data, bits);
        }//end blocks (while loop)

        //mean extrinsic mutual info over all blocks
        extrinsic_mutual_info(en) /= nb_blocks_lim;
        extrinsic_mutual_info_p(en) /= nb_blocks_lim;
    }

    cout << "IA = " << fixed << round_to_infty(apriori_mutual_info) << endl;
    cout << "IE = " << round_to_infty(extrinsic_mutual_info) << endl;
    cout << "IE^p = " << round_to_infty(extrinsic_mutual_info_p) << endl;


}

IGNORE_TEST(HelloItpp,Fileter)
{

  /*-----------------------------------------------------------------------------
   *  test of filter
   *-----------------------------------------------------------------------------*/
  cout << "======================================" << endl;
  cout << "        Test of filter routines       " << endl;
  cout << "======================================" << endl;

  // Test signals
  vec x = randn(20), x2 = randn(20);
  cvec cx = randn_c(20), cx2 = randn_c(20);

  cout << "Input signals: " << endl;
  cout << "x = " << x << endl;
  cout << "x2 = " << x2 << endl;
  cout << "cx = " << cx << endl;
  cout << "cx2 = " << cx2 << endl;

  // Filter coefficients
  vec b(10);
  b.ones();
  b(0) += 0.1;
  cvec cb(10);
  cb.ones();
  cb(0) += 0.1;

  cvec ca(2);
  ca(0) = 1.0;
  ca(1) = -1.0;
  vec a(2);
  a(0) = 1.0;
  a(1) = -1.0;

  cout << "Filter coefficients: " << endl;
  cout << "b = " << b << endl;
  cout << "cb = " << cb << endl;
  cout << "a = " << a << endl;
  cout << "ca = " << ca << endl;

  vec y, y2;
  vec s1, s2;

  cvec cy, cy2;
  cvec cs1, cs2;

  cout << endl << "-----------------------------------------------" << endl;
  cout << "MA Filter: " << endl;

  MA_Filter<double, double, double> H(b);
  MA_Filter<complex<double>, double, complex<double> > CH(b);
  MA_Filter<complex<double>, complex<double>, complex<double> > C(cb);

  y = H(x);
  s1 = H.get_state();
  cout << "y = " << y << endl;
  cout << "s1 = " << s1 << endl;
  y2 = H(x2);
  s2 = H.get_state();
  cout << "y2 = " << y2 << endl;
  cout << "s2 = " << s2 << endl;
  cout << "Redo, reseting to state s1" << endl;
  H.set_state(s1);
  y2 = H(x2);
  s2 = H.get_state();
  cout << "y2 = " << y2 << endl;
  cout << "s2 = " << s2 << endl;

  cy = CH(cx);
  cs1 = CH.get_state();
  cout << "cy = " << cy << endl;
  cout << "cs1 = " << cs1 << endl;
  cy2 = CH(cx2);
  cs2 = CH.get_state();
  cout << "cy2 = " << cy2 << endl;
  cout << "cs2 = " << cs2 << endl;
  cout << "Redo, reseting to state s1" << endl;
  CH.set_state(cs1);
  cy2 = CH(cx2);
  cs2 = CH.get_state();
  cout << "cy2 = " << cy2 << endl;
  cout << "cs2 = " << cs2 << endl;

  cy = C(cx);
  cout << "cy = " << cy << endl;

  y =  filter(b, 1, x);
  cout << "y = " << y << endl;

  cy =  filter(b, 1, cx);
  cout << "cy = " << cy << endl;

  cy =  filter(cb, 1, cx);
  cout << "cy = " << cy << endl;


  cout << endl << "-----------------------------------------------" << endl;
  cout << "AR Filter: " << endl;

  AR_Filter<double, double, double> HAR(a);
  AR_Filter<complex<double>, double, complex<double> > CHAR(a);
  AR_Filter<complex<double>, complex<double>, complex<double> > CAR(ca);

  y = HAR(x);
  s1 = HAR.get_state();
  cout << "y = " << y << endl;
  cout << "s1 = " << s1 << endl;
  y2 = HAR(x2);
  s2 = HAR.get_state();
  cout << "y2 = " << y2 << endl;
  cout << "s2 = " << s2 << endl;
  cout << "Redo, reseting to state s1" << endl;
  HAR.set_state(s1);
  y2 = HAR(x2);
  s2 = HAR.get_state();
  cout << "y2 = " << y2 << endl;
  cout << "s2 = " << s2 << endl;

  cy = CHAR(cx);
  cs1 = CHAR.get_state();
  cout << "cy = " << cy << endl;
  cout << "cs1 = " << cs1 << endl;
  cy2 = CHAR(cx2);
  cs2 = CHAR.get_state();
  cout << "cy2 = " << cy2 << endl;
  cout << "cs2 = " << cs2 << endl;
  cout << "Redo, reseting to state s1" << endl;
  CHAR.set_state(cs1);
  cy2 = CHAR(cx2);
  cs2 = CHAR.get_state();
  cout << "cy2 = " << cy2 << endl;
  cout << "cs2 = " << cs2 << endl;

  cy = CAR(cx);
  cout << "cy = " << cy << endl;

  y =  filter(1, a, x);
  cout << "y = " << y << endl;

  cy =  filter(1, a, cx);
  cout << "cy = " << cy << endl;

  cy =  filter(1, ca, cx);
  cout << "cy = " << cy << endl;


  cout << endl << "-----------------------------------------------" << endl;
  cout << "ARMA Filter: " << endl;

  ARMA_Filter<double, double, double> HARMA(b, a);
  ARMA_Filter<complex<double>, double, complex<double> > CHARMA(b, a);
  ARMA_Filter<complex<double>, complex<double>, complex<double> > CARMA(cb, ca);

  y = HARMA(x);
  s1 = HARMA.get_state();
  cout << "y = " << y << endl;
  cout << "s1 = " << s1 << endl;
  y2 = HARMA(x2);
  s2 = HARMA.get_state();
  cout << "y2 = " << y2 << endl;
  cout << "s2 = " << s2 << endl;
  cout << "Redo, reseting to state s1" << endl;
  HARMA.set_state(s1);
  y2 = HARMA(x2);
  s2 = HARMA.get_state();
  cout << "y2 = " << y2 << endl;
  cout << "s2 = " << s2 << endl;

  cy = CHARMA(cx);
  cs1 = CHARMA.get_state();
  cout << "cy = " << cy << endl;
  cout << "cs1 = " << cs1 << endl;
  cy2 = CHARMA(cx2);
  cs2 = CHARMA.get_state();
  cout << "cy2 = " << cy2 << endl;
  cout << "cs2 = " << cs2 << endl;
  cout << "Redo, reseting to state s1" << endl;
  CHARMA.set_state(cs1);
  cy2 = CHARMA(cx2);
  cs2 = CHARMA.get_state();
  cout << "cy2 = " << cy2 << endl;
  cout << "cs2 = " << cs2 << endl;

  cy = CARMA(cx);
  cout << "cy = " << cy << endl;

  y =  filter(b, a, x);
  cout << "y = " << y << endl;

  cy =  filter(b, a, cx);
  cout << "cy = " << cy << endl;

  cy =  filter(cb, ca, cx);
  cout << "cy = " << cy << endl;

}


IGNORE_TEST(HelloItpp,BlerBer)
{
  cout << "==========================" << endl;
  cout << "    BERC and BLER test    " << endl;
  cout << "==========================" << endl << endl;

  const int block_size = 20;
  const ivec error_pos = "3 5 9 17";
  BERC berc;
  BLERC blerc(block_size);

  for (int i = 0; i < 100; ++i) {
    bvec input = randb(block_size);
    bvec output = input;

    // introduce some errors
    if (i < 80)
      for (int j = 0; j < error_pos.size(); ++j)
        output(error_pos(j)) = !output(error_pos(j));
    // extend the output vector by one bit
    output = concat(output, bin(1));
    // count errors
    berc.count(input, output);
    blerc.count(input, output);
  }
  cout << "BLER = " << blerc.get_errorrate() << endl
       << "Block errors = " << blerc.get_errors() << endl
       << "Correct blocks = " << blerc.get_corrects() << endl
       << "Total blocks = " << blerc.get_total_blocks() << endl << endl;
  blerc.clear();

  cout << "BER = " << berc.get_errorrate() << endl
       << "Bit errors = " << berc.get_errors() << endl
       << "Correct bits = " << berc.get_corrects() << endl
       << "Total bits = " << berc.get_total_bits() << endl;


}

IGNORE_TEST(HelloItpp,Eigen)
{
  cout << "================================" << endl;
  cout << "  Test of eigenvalue routines   " << endl;
  cout << "================================" << endl;

  const double thresh = 1e-13;

  {
    cout << "Real symmetric matrix" << endl;
    mat A = randn(5, 5);
    A = transpose(A) * A; // make it symmetic
    mat V;
    vec d;
    eig_sym(A, d, V);

    cout << "A = " << round_to_zero(A) << endl;
    cout << "norm(A*V-V*diag(d)) = "
         << round_to_zero(norm(A * V - V * diag(d)), thresh) << endl;
  }

  {
    cout << endl << "Real non-symmetric matrix" << endl;
    mat A = randn(5, 5);
    cmat V;
    cvec d;
    eig(A, d, V);

    cout << "A = " << round_to_zero(A) << endl;
    cout << "norm(A*V-V*diag(d)) = "
         << round_to_zero(norm(A * V - V * diag(d)), thresh) << endl;
  }

  {
    cout << endl << "Complex hermitian matrix" << endl;
    cmat A = randn_c(5, 5);
    A = transpose(conj(A)) * A; // make it hermitian
    cmat V;
    vec d;
    eig_sym(A, d, V);

    cout << "A = " << round_to_zero(A) << endl;
    cout << "norm(A*V-V*diag(d)) = "
         << round_to_zero(norm(A * V - V * diag(d)), thresh) << endl;
  }

  {
    cout << endl << "Complex non-hermitian matrix" << endl;
    cmat A = randn_c(5, 5);
    cmat V;
    cvec d;
    eig(A, d, V);

    cout << "A = " << round_to_zero(A) << endl;
    cout << "norm(A*V-V*diag(d)) = "
         << round_to_zero(norm(A * V - V * diag(d)), thresh) << endl;
  }


}

IGNORE_TEST(HelloItpp, Fastica)
{
  FILE * fpin = NULL;
  float tmp = 0.0;

  // Separate nrIC independent components in nrSamples samples
  int nrSamples = 0, nrIC = 0;

  fpin = fopen(FASTICA_TEST_FILE, "r");
  if (fpin == 0) {
    cerr << "Error: Could not open FASTICA_TEST_FILE for reading" << endl;
    return ;
  }

  cout << "=====================================" << endl;
  cout << "   Test program for FastICA / IT++   " << endl;
  cout << "=====================================" << endl;

  int ret = fscanf(fpin, "%d", &nrSamples);
  ret = fscanf(fpin, "%d", &nrIC);
  cout <<"ret of fscan FastICA" << ret << endl;
  mat X = zeros(nrIC, nrSamples);

  for (int i = 0; i < nrSamples; i++)
    for (int j = 0; j < nrIC; j++) {
      ret = fscanf(fpin , "%f", &tmp);
      X(j, i) = tmp;
    }

  fclose(fpin);

  // Instantiate an ICA object with default parameters : SYMM approach and
  // POW3 non-linearity
  // Be sure that :
  // - nrSamples = number of samples = nb of columns of the input matrix
  // - nrIC = number of sensors = nb of rows of the input matrix
  Fast_ICA my_fastica(X);

  // Set number of independent components to separate :
  // By default, this value is taken from the dimension of
  // the input data. This line is for illustration purposes.
  // May help in some cases.
  my_fastica.set_nrof_independent_components(nrIC);

  // Perform ICA
  my_fastica.separate();

  cout << "Use default parameters:" << endl;

  // Get results : mixing and separating matrices
  cout << "Mixing matrix = " << my_fastica.get_mixing_matrix() << endl;
  cout << "Separation matrix = " << my_fastica.get_separating_matrix() << endl;

  // Get result : separated independent components
  cout << endl << "separated independent components = "
       << my_fastica.get_independent_components();

  // Another test with other parameters
  cout << "==========================================================" << endl;
  cout << "Use Gaussian non-linearity and deflation approach :" << endl;

  Fast_ICA my_fastica2(X);

  // Set GAUSS non-linearity
  my_fastica2.set_non_linearity(FICA_NONLIN_GAUSS);

  // Use deflation approach : IC are computed one by one
  my_fastica2.set_approach(FICA_APPROACH_DEFL);

  // Perform ICA
  my_fastica2.separate();

  // Get results
  cout << "Mixing matrix = " << my_fastica2.get_mixing_matrix() << endl;
  cout << "Separation matrix = " << my_fastica2.get_separating_matrix() << endl;
  cout << endl << "separated independent components = "
       << my_fastica2.get_independent_components() << endl;

  cout << "End of Fast_ICA test execution. " << endl;


}

IGNORE_TEST(HelloItpp, Fastmath)
{
  cout << "====================================================" << endl;
  cout << "              Test of fastmath" << endl;
  cout << "====================================================" << endl;

  mat m0("1 2 3;4 5 6;7 8 9"), mv0("2;3;1");
  vec v0("2 3 1");

  cout << "sub_v_vT_m: the slow and fast way" << endl;
  cout << (m0 - mv0*transpose(mv0)*m0) << endl;
  sub_v_vT_m(m0, v0);
  cout << m0 << endl;

  m0 = "1 2 3;4 5 6;7 8 9";

  cout << endl << "sub_m_v_vT: the slow and fast way" << endl;
  cout << (m0 - m0*mv0*transpose(mv0)) << endl;
  sub_m_v_vT(m0, v0);
  cout << m0 << endl;


}

IGNORE_TEST(HelloItpp, Determine)
{
  cout << "===================================" << endl;
  cout << "    Test of Determinant routines   " << endl;
  cout << "===================================" << endl;

  {
    cout << "Real matrix" << endl;
    mat X = randn(5, 5);
    double d;
    d = det(X);
    cout << "X = " << round_to_zero(X) << endl;
    cout << "det(X) = " << round_to_zero(d) << endl;

    X = randn(5, 5);
    d = det(X);
    cout << "X = " << round_to_zero(X) << endl;
    cout << "det(X) = " << round_to_zero(d) << endl;
  }

  {
    cout << endl << "Complex matrix" << endl;
    cmat X = randn_c(5, 5);
    complex<double> d;
    d = det(X);
    cout << "X = " << round_to_zero(X) << endl;
    cout << "det(X) = " << round_to_zero(d) << endl;

    X = randn_c(5, 5);
    d = det(X);
    cout << "X = " << round_to_zero(X) << endl;
    cout << "det(X) = " << round_to_zero(d) << endl;
  }


}

IGNORE_TEST(HelloItpp, Convert)
{
  cout << "==================================" << endl
       << "   Test of conversion functions" << endl
       << "==================================" << endl;

  vec v = "-4:0.5:5.5";
  v = concat(v, randn(10) * 20 - 10);

  cout << "v = " << v << endl << endl;

  cout << "round(v) = " << round_to_zero(round(v)) << endl;
  cout << "round_i(v) = " << round_i(v) << endl << endl;

  cout << "ceil(v) = " << round_to_zero(ceil(v)) << endl;
  cout << "ceil_i(v) = " << ceil_i(v) << endl << endl;

  cout << "floor(v) = " << round_to_zero(floor(v)) << endl;
  cout << "floor_i(v) = " << floor_i(v) << endl << endl;

  bvec b = randb(15);

  cout << "b = " << b << endl << endl;

  int i = bin2dec(b);
  cout << "bin2dec(b) = " << i << endl;
  cout << "dec2bin(bin2dec(b)) = " << dec2bin(i) << endl;
  i = bin2dec(b, false);
  cout << "bin2dec(b, false) = " << i << endl;
  cout << "dec2bin(bin2dec(b, false), false) = " << dec2bin(i, false)
       << endl << endl;

  ivec iv = bin2oct(b);
  cout << "bin2oct(b) = " << iv << endl;
  cout << "oct2bin(bin2oct(b)) = " << oct2bin(iv) << endl;
  cout << "oct2bin(bin2oct(b), 1) = " << oct2bin(iv, 1) << endl << endl;

  iv = bin2pol(b);
  cout << "bin2pol(b) = " << iv << endl;
  cout << "pol2bin(bin2pol(b)) = " << pol2bin(iv) << endl << endl;


}
IGNORE_TEST(HelloItpp,Conv)
{
  cout << "====================================" << endl;
  cout << "    Test of convolutional coders    " << endl;
  cout << "====================================" << endl;

  Array<ivec> spectrum, spectrum_fast, spectrum_punct, spectrum_punct_fast;
  spectrum.set_size(2);
  spectrum_fast.set_size(2);
  spectrum_punct.set_size(2);
  spectrum_punct_fast.set_size(2);

  Convolutional_Code code;
  Punctured_Convolutional_Code code_punct;
  BPSK bpsk;
  BERC berc;

  const int no_bits = 2500;
  const int packet_size = 500;

  int coded_packet_size;
  bvec bits, tail_coded_bits, tail_decoded_bits, tailbite_coded_bits,
  tailbite_decoded_bits, trunc_coded_bits, trunc_decoded_bits;
  vec symbols;
  ivec dist_profile;

  ivec G(2);
  G(0) = 0133;
  G(1) = 0171;
  int L = max(int2bits(G(0)), int2bits(G(1))); // L = 7

  code.set_generator_polynomials(G, L);
  code_punct.set_generator_polynomials(G, L);

  bmat punct_matrix = "1 0 1; 1 1 0"; // results in R = 3/4
  code_punct.set_puncture_matrix(punct_matrix);


  cout << "------------------------------------------------------------------------------" << endl;
  cout << "1) Rate 1/2 code" << endl;
  cout << "------------------------------------------------------------------------------" << endl;

  cout << "Catastrophic test = " << code.catastrophic() << endl;
  cout << "Code rate         = " << code.get_rate() << endl << endl;

  code.calculate_spectrum(spectrum, 10, 10);
  code.fast(spectrum_fast, 10, 10);
  code.distance_profile(dist_profile, 10);

  cout << "Spectrum:" << endl;
  cout << "* Ad = " << spectrum(0) << endl;
  cout << "* Cd = " << spectrum(1) << endl;

  cout << "Spectrum, fast:" << endl;
  cout << "* Ad = " << spectrum_fast(0) << endl;
  cout << "* Cd = " << spectrum_fast(1) << endl << endl;

  cout << "Distance profile  = " << dist_profile << endl << endl;

  cout << "Tail method test. Printing 30 bits starting from bit 1400:" << endl;
  bits = randb(no_bits);
  cout << "* Input bits    = " << bits.mid(1400, 30) << endl;
  tail_coded_bits = code.encode_tail(bits);
  cout << "* Coded bits    = " << tail_coded_bits.mid(1400, 30) << endl;
  bpsk.modulate_bits(tail_coded_bits, symbols);
  tail_decoded_bits = code.decode_tail(symbols);
  cout << "* Decoded bits  = " << tail_decoded_bits.mid(1400, 30) << endl;
  berc.count(bits, tail_decoded_bits);
  cout << "BER = " << berc.get_errorrate() << endl << endl;

  cout << "Tailbite method test. Printing 30 bits starting from bit 1400:"
       << endl;
  cout << "* Input bits    = " << bits.mid(1400, 30) << endl;
  tailbite_coded_bits = code.encode_tailbite(bits);
  cout << "* Coded bits    = " << tailbite_coded_bits.mid(1400, 30) << endl;
  bpsk.modulate_bits(tailbite_coded_bits, symbols);
  tailbite_decoded_bits = code.decode_tailbite(symbols);
  cout << "* Decoded bits  = " << tailbite_decoded_bits.mid(1400, 30) << endl;
  berc.clear();
  berc.count(bits, tailbite_decoded_bits);
  cout << "BER = " << berc.get_errorrate() << endl << endl;

  cout << "Trunc method test. Printing 30 bits starting from bit 1400:" << endl;
  cout << "* Input bits    = " << bits.mid(1400, 30) << endl;
  trunc_coded_bits.set_size(0);
  for (int i = 0; i < no_bits / packet_size; i++) {
    trunc_coded_bits = concat(trunc_coded_bits,
                              code.encode_trunc(bits.mid(i * packet_size,
                                                         packet_size)));
  }
  cout << "* Coded bits    = " << trunc_coded_bits.mid(1400, 30) << endl;
  bpsk.modulate_bits(trunc_coded_bits, symbols);
  trunc_decoded_bits.set_size(0);
  coded_packet_size = round_i(packet_size / code.get_rate());
  for (int i = 0; i < no_bits / packet_size; i++) {
    trunc_decoded_bits =
      concat(trunc_decoded_bits,
             code.decode_trunc(symbols.mid(i * coded_packet_size,
                                           coded_packet_size)));
  }
  cout << "* Decoded bits  = " << trunc_decoded_bits.mid(1400, 30) << endl;
  berc.clear();
  berc.count(bits, trunc_decoded_bits);
  cout << "BER = " << berc.get_errorrate() << endl << endl;


  cout << "------------------------------------------------------------------------------" << endl;
  cout << "1) Punctured code (R = 3/4)" << endl;
  cout << "------------------------------------------------------------------------------" << endl;

  cout << "Catastrophic test = " << code_punct.catastrophic() << endl;
  cout << "Code rate         = " << code_punct.get_rate() << endl;
  cout << "Puncture matrix   = " << code_punct.get_puncture_matrix() << endl
       << endl;

  cout << "Tail method test. Printing 30 bits starting from bit 1400:" << endl;
  bits = randb(no_bits);
  cout << "* Input bits    = " << bits.mid(1400, 30) << endl;
  tail_coded_bits = code_punct.encode_tail(bits);
  cout << "* Coded bits    = " << tail_coded_bits.mid(1400, 30) << endl;
  bpsk.modulate_bits(tail_coded_bits, symbols);
  tail_decoded_bits = code_punct.decode_tail(symbols);
  cout << "* Decoded bits  = " << tail_decoded_bits.mid(1400, 30) << endl;
  berc.count(bits, tail_decoded_bits);
  cout << "BER = " << berc.get_errorrate() << endl << endl;

  cout << "Tailbite method test. Printing 30 bits starting from bit 1400:"
       << endl;
  cout << "* Input bits    = " << bits.mid(1400, 30) << endl;
  tailbite_coded_bits = code_punct.encode_tailbite(bits);
  cout << "* Coded bits    = " << tailbite_coded_bits.mid(1400, 30) << endl;
  bpsk.modulate_bits(tailbite_coded_bits, symbols);
  tailbite_decoded_bits = code_punct.decode_tailbite(symbols);
  cout << "* Decoded bits  = " << tailbite_decoded_bits.mid(1400, 30) << endl;
  berc.clear();
  berc.count(bits, tailbite_decoded_bits);
  cout << "BER = " << berc.get_errorrate() << endl << endl;

  cout << "Trunc method test. Printing 30 bits starting from bit 1400:" << endl;
  cout << "* Input bits    = " << bits.mid(1400, 30) << endl;
  trunc_coded_bits.set_size(0);
  for (int i = 0; i < no_bits / packet_size; i++) {
    trunc_coded_bits = concat(trunc_coded_bits,
                              code_punct.encode_trunc(bits.mid(i * packet_size,
                                                               packet_size)));
  }
  cout << "* Coded bits    = " << trunc_coded_bits.mid(1400, 30) << endl;
  bpsk.modulate_bits(trunc_coded_bits, symbols);
  trunc_decoded_bits.set_size(0);
  coded_packet_size = round_i(packet_size / code_punct.get_rate());
  for (int i = 0; i < no_bits / packet_size; i++) {
    trunc_decoded_bits =
      concat(trunc_decoded_bits,
             code_punct.decode_trunc(symbols.mid(i * coded_packet_size,
                                                 coded_packet_size)));
  }
  cout << "* Decoded bits  = " << trunc_decoded_bits.mid(1400, 30) << endl;
  berc.clear();
  berc.count(bits, trunc_decoded_bits);
  cout << "BER = " << berc.get_errorrate() << endl << endl;


}

IGNORE_TEST(HelloItpp,Hamming)
{
	cout << graycode(5) << endl;

	bvec b1 = randb(10);
	bvec b2 = randb(10);
	cout << hamming_distance(b1, b2) << endl;

	cout << weight(b1) << endl;

	cout << waterfilling("1 2 3", 0.001) << endl;
	cout << waterfilling("1 2 3", 0.01) << endl;
	cout << waterfilling("1 2 3", 0.1) << endl;
	cout << waterfilling("1 2 3", 0.5) << endl;
	cout << waterfilling("1 2 3", 1) << endl;
	cout << waterfilling("1 2 3", 5) << endl;
	cout << waterfilling("1 2 3", 100) << endl;

}

IGNORE_TEST(HelloItpp,CircularBuffer)
{
	cout << "\nHelloItpp CircularBuffer\n" << endl;
	//Set a fixed seed to the random number generator
	RNG_reset(12345);

	int nrof_elements;
	ivec index_vec;
	vec a = randn(3);
	vec b = randn(8);
	vec out_vec;
	double out;
	Array <double> out_array;

	Circular_Buffer<double> cb1(10);

	//Put the elements of a to the buffer
	cb1.put(a);

	//Vector output: Peek at the two oldest elements of the buffer (the two first extracted)
	cb1.peek(out_vec, 2);
	cout << "peek(out_vec,2) = " << out_vec << ": display 2 first elements of the buffer, without affecting the content" << endl;

	//Vector output: Peek at all elements of the buffer in reverse order
	cb1.peek_reverse(out_vec);
	cout << "peek_reverse(out_vec,-1) = " << out_vec << ": display buffer, without affecting the content" << endl;

	//Array output: Peek at all elements of the buffer in reverse order
	cb1.peek_reverse(out_array);
	cout << "peek_reverse(out_array,-1) = " << out_array << ": display buffer, without affecting the content" << endl;

	//Put the elements of \c b to the buffer
	cb1.put(b);

	//Extract the oldest element of the buffer
	cb1.get(out_vec, 1);
	cout << "get(out_vec,1) = " << out_vec << endl ;

	//Extract all element of the buffer
	cb1.get(out_vec);
	cout << "get(out_vec) = " << out_vec << endl ;

	cb1.get(out_vec);
	cb1.get(out_array, 0);
	cout << "get(out_vec) = " << out_vec;
	cout << "; get(out_array,0) = " << out_array << ": empty buffer, no content" << endl;

	for (int i = 0;i < a.length();i++) { cb1.put(a(i)); }
	for (int i = 0;i < b.length();i++) { cb1.put(b(i)); }

	cout << "buffer size = " << cb1.size() << endl;
	cout << "nrof_elements = " << cb1.nrof_elements() << endl;

	nrof_elements = cb1.nrof_elements();
	for (int i = 0;i < nrof_elements;i++) {
		cout << "i = " << i;
		cout << ": peek() = " << cb1.peek();
		cout << "; peek_reverse() = " << cb1.peek_reverse();
		cout << "; get() = " << cb1.get();
		cout << endl;
	}

	//Test of error handling
	//cout << "get() = " << cb1.get() << endl;

	cb1.put(b);
	cout << "buffer size = " << cb1.size() << endl;
	cout << "nrof_elements = " << cb1.nrof_elements() << endl;

	cb1.peek(out_vec, -1);
	cout << "peek(out_vec) = " << out_vec << endl;

	index_vec = "5 3 7 1";
	for (int i = 0;i < index_vec.size();i++) {
		cb1.peek(index_vec(i), out);
		cout << "peek at index " << index_vec(i);
		cout << ": peek(" << index_vec(i) << ") = " << out << endl;
	}

	cb1.peek(index_vec, out_vec);
	cout << "peek at indices " << index_vec << ":" << endl;
	cout << "peek(index_vec,out_vec) = " << out_vec << endl;

//	cb1.set_size(15, true);
//	cout << "buffer size = " << cb1.size() << endl;
//	cout << "nrof_elements = " << cb1.nrof_elements() << endl;
//
//	cb1.peek(out_vec, -1);
//	cout << "peek(out_vec) = " << out_vec << endl;

//	cb1.set_size(5, true);
//	cout << "buffer size = " << cb1.size() << endl;
//	cout << "nrof_elements = " << cb1.nrof_elements() << endl;
//
//	cb1.peek(out_vec, -1);
//	cout << "peek(out_vec) = " << out_vec << endl;
//
	cout << "Copy constructor: " << endl;
	Circular_Buffer<double> cb2(cb1);
	cb2.peek(out_vec, -1);
	cout << "Circular_Buffer<double> cb2(cb1): cb2.peek(out_vec) = " << out_vec << endl;

	cout << "Copy operator: " << endl;
	Circular_Buffer<double> cb3 = cb1;
	cb3.peek(out_array, -1);
	cout << "Circular_Buffer<double> cb3=cb1: cb3.peek(out_array) = " << out_array << endl;

}

IGNORE_TEST(HelloItpp,Cholesky)
{
	cout << "================================" << endl;
	cout << " Test of Cholesky routines " << endl;
	cout << "================================" << endl;

	{
	cout << "Real matrix" << endl;
	mat X, F;
	bool ok = false;

		while (!ok) {
		X = randn(5, 5);
		X = transpose(X) * X; // create a symmetric matrix
		// Make diagonal real and positive
		for (int i = 0; i < X.cols(); i++)
		X(i, i) = std::abs(X(i, i));

		ok = chol(X, F);
		cout << "X = " << round_to_zero(X) << endl;
		if (!ok)
			cout << "matrix is not positive definite" << endl;
		else {
			cout << "norm(X - F^T*F) = "
			<< round_to_zero(norm(X - transpose(F) * F)) << endl;
		}
	}
	}

	{
	cout << endl << "Complex matrix" << endl;
	cmat X, F;
	bool ok = false;

	while (!ok) {
		X = randn_c(5, 5);
		X = hermitian_transpose(X) * X; // create a symmetric matrix
		// Make diagonal real and positive
		for (int i = 0; i < X.cols(); i++)
			X(i, i) = std::abs(real(X(i, i)));

		ok = chol(X, F);
		cout << "X = " << round_to_zero(X) << endl;

		if (!ok)
			cout << "matrix is not positive definite" << endl;
		else {
			cout << "norm(X - F^H*F) = "
			<< round_to_zero(norm(X - hermitian_transpose(F) * F)) << endl;
		}
	}
}

}

IGNORE_TEST(HelloItpp,Blas)
{
	cout << "\n HelloItpp,Blas\n" << endl;


	cout.setf(ios::fixed);
	cout.precision(4);

	// dot() tests
	{
	vec a = randn(10);
	vec b = randn(10);
	cout << "a = " << a << endl;
	cout << "b = " << b << endl;
	cout << "dot(a,b) = " << dot(a, b) << endl << endl;

	cvec c = randn_c(6);
	cvec d = randn_c(6);
	cout << "c = " << c << endl;
	cout << "d = " << d << endl;
	cout << "dot(c,d) = " << dot(c, d) << endl << endl;
	}

	// outer_product() tests
	{
	vec a = randn(4);
	vec b = randn(7);
	cout << "a = " << a << endl;
	cout << "b = " << b << endl;
	cout << "outer_product(a,b) = " << outer_product(a, b) << endl;

	cvec c = randn_c(4);
	cvec d = randn_c(7);
	cout << "c = " << c << endl;
	cout << "d = " << d << endl;
	cout << "outer_product(c,d) = " << outer_product(c, d) << endl;
	cout << "outer_product(c,d,true) = " << outer_product(c, d, true) << endl << endl;
	}

	// Mat *= Mat operator test
	{
	mat M1 = randn(3, 5);
	mat N1 = randn(5, 2);
	cout << "M = " << M1 << endl;
	cout << "N = " << N1 << endl;
	M1 *= N1;
	cout << "M *= N;\nM = " << M1 << endl << endl;

	cmat M2 = randn_c(4, 4);
	cmat N2 = randn_c(4, 2);
	cout << "M = " << M2 << endl;
	cout << "N = " << N2 << endl;
	M2 *= N2;
	cout << "M *= N;\nM = " << M2 << endl << endl;
	}

	// Vec = Mat * Vec operator test
	{
	mat M1 = randn(3, 4);
	vec v1 = randn(4);
	cout << "M = " << M1 << endl;
	cout << "v = " << v1 << endl;
	cout << "out = M * v = " << M1 * v1 << endl << endl;

	cmat M2 = randn_c(3, 2);
	cvec v2 = randn_c(2);
	cout << "M = " << M2 << endl;
	cout << "v = " << v2 << endl;
	cout << "out = M * v = " << M2 * v2 << endl << endl;
	}

}

IGNORE_TEST(HelloItpp, Channel)
{
	cout.setf(ios::fixed);
	cout.precision(3);

	double Ts = 100e-9; // channel sampling time Ts = 100ns
	double fd_norm = 100e3 * Ts; // normalised Doppler fd_norm = 10kHz * Ts

	Channel_Specification cost207ra(COST207_RA);
	Channel_Specification cost207bu(COST207_BU);

	TDL_Channel tdl_207ra(cost207ra, Ts);
	TDL_Channel tdl_207bu(cost207bu, Ts);

	cmat ch_coeffs;

	tdl_207ra.set_fading_type(Independent);
	tdl_207ra.generate(10, ch_coeffs);
	cout << "Independent Fading Generator\n"
	   << "----------------------------\n"
	   << ch_coeffs << "\n\n";

	tdl_207bu.set_fading_type(Static);
	tdl_207bu.generate(5, ch_coeffs);
	cout << "Static Fading Generator\n"
	   << "-----------------------\n"
	   << ch_coeffs << "\n\n";

	tdl_207bu.set_fading_type(Correlated);
	tdl_207bu.set_norm_doppler(fd_norm);
	tdl_207bu.generate(10, ch_coeffs);
	cout << "Correlated Fading Generator (Rice method)\n"
	   << "-----------------------------------------\n"
	   << round_to_zero(ch_coeffs) << "\n\n";

	tdl_207ra.set_correlated_method(FIR);
	tdl_207ra.set_norm_doppler(fd_norm);
	tdl_207ra.generate(10, ch_coeffs);
	cout << "Correlated Fading Generator (FIR method)\n"
	   << "----------------------------------------\n"
	   << round_to_zero(ch_coeffs) << "\n\n";

	tdl_207ra.set_correlated_method(IFFT);
	tdl_207ra.generate(200, ch_coeffs);
	cout << "Correlated Fading Generator (IFFT method)\n"
	   << "-----------------------------------------\n"
	   << round_to_zero(ch_coeffs.get_rows(0, 9)) << "\n\n";

	return ;
//	STRCMP_EQUAL("Hello World!\n", buffer->asCharString());
}
